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We calculate how strongly one can put constraints on alternative theories of gravity such as 
Brans-Dicke and massive graviton theories with LISA. We consider inspiral gravitational waves from 
a compact binary composed of a neutron star (NS) and an intermediate mass black hole (IMBH) 
in Brans-Dicke (BD) theory and that composed of 2 super massive black holes (SMBHs) in massive 
graviton theories. We use the restricted 2PN waveforms including the effects of spins. We also take 
, both precession and eccentricity of the orbit into account. For simplicity, we set the fiducial value for 

. the spin of one of the binary constituents to zero so that we can apply the approximation called simple 

precession. We perform the Monte Carlo simulations of 10* binaries, estimating the determination 
accuracy of binary parameters including the BD parameter ujbi) and the Compton wavelength of 
graviton Ag for each binary using the Fisher matrix method. We find that including both the spin- 
spin coupling a and the eccentricity e into the binary parameters reduces the determination accuracy 
by an order of magnitude for the Brans-Dicke case, whilst it has less influence on massive graviton 
theories. On the other hand, including precession enhances the constraint on ujbd only 20% but it 
increases the constraint on Xg by an order of magnitude. Using a (1.4 + 1000)Mq NS/BH binary 
of SNR=\/200, one can put a constraint a;BD > 6944, whilst using a (10^ + 1O'^)M0 BH/BH binary 
O at 3Gpc, one can put Ag > 3.06 x lO^^cm, on average. The latter is 4 orders of magnitude stronger 

than the one obtained from the solar system experiment. These results are consistent with previous 
»' I ' results within uncontrolled errors and indicate that the effects of precession and eccentricity must 

{3J[), be taken carefully in the parameter estimation analysis. 

PACS numbers: Valid PACS appear here 

> 

(N 

■"si" ' Recent observations (e.g. type la SNe [1]) show that the expansion of the Universe is accelerating at the mo- 

\^ : ment. One possible way to explain this accelerated expansion is to introduce some undiscovered matter field such 
as quintessence or k-essence (for a recent review, see [2]). Another possibility is to modify gravitational theory from 
0^ . general relativity. The easiest modification is to add some scalar degrees of freedom to gravity. This type of theory 
■ is called the scalar-tensor theory [3]. This theory also appears to be a candidate for solving the inflation problem, 
J>. I known as the hyper-extended inflation [4]. Furthermore, scalar fields called dilatons and moduli may play the role 
■ '~j ■ of these scalar degrees of freedom in the context of the string theory. Therefore it is very important to put a strong 
rN observational constraint on this theory. In this paper, we focus on Brans-Dicke (BD) theory [5] as a representative of 
^ I scalar-tensor theory. This theory has the so-called Brans-Dicke parameter wbd and by taking the limit wbd — >■ oo, it 
■ " " ' reduces to general relativity. The current strongest constraint on the Brans-Dicke parameter is wbd > 40000, which 
is obtained from the solar system experiment, measuring the Shapiro time delay using the Saturn probe satellite 
Cassini [6]. 

Another type of modified theory of gravity introduces a finite mass mg to a graviton (see [7] for a recent review). 
This type of theory is called the massive graviton theory. Originally, Fierz and Pauli [8] proposed a Lorentz-invariant 
massive gravity by simply adding a graviton mass term to the Einstein-Hilbert action at the quadratic level. However, 
it was found that the linearised Fierz-Pauli theory does not approach linearised general relativity in the massless limit. 
This is the so-called van Dam-Veltman-Zakharov (vDVZ) discontinuity [9, 10]. This is due to the fact that the helicity- 
component of the graviton does not decouple from matter. This discontinuity seems to contradict with solar system 
experiments, but Vainshtein pointed out that the effect of nonlinearity might be important [11]. In general relativity, 
the linear approximation is valid for distances much larger than the Schwarzschild radius of the source. However, in 
massive gravity, linear approximation breaks down already at a distance much longer than the Schwarzschild radius. 
Indeed, this Vainshtein effect has been shown to work in the DGP braneworld model [12, 13]. Also, Rubakov [14] 
and Dubovsky [15] have proposed Lorentz violating massive gravity theories which evade pathologies related to the 
vDVZ discontinuity. Therefore it is also important to put constraint on the graviton mass (or the graviton Compton 
wavelength Xg = h/mgc). When the graviton is massive, the gravitational potential is modified from Newtonian to 
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Yukawa type. Then, the effective gravitational constant depends on the distance from the gravitational source, which 
changes the Kepler's third law from Newtonian. Verification of this third law in the solar system experiment puts a 
lower bound on the Compton wavelength as Xg > 2.8 x lO^^cm [16]. 

In this paper, we estimate the possible constraint we can get on wbd and by detecting gravitational waves from the 
inspiral of processing eccentric compact binaries with LISA [17, 18]. The constraint on wbd using binary gravitational 
waves was first discussed by Eardley [19]. In Brans-Dicke theory, the additional scalar field contains the dipole 
radiation [20, 21]. This modifies the binary's orbital evolution from the one in general relativity. Eardley mentioned 
that it is possible to constrain wbd by measuring the rate of the secular decrease in the orbital period. The change in 
the orbital evolution due to this dipole radiation modifies the phasing of the gravitational waveform. Will [22] applied 
the matched filtering analysis and calculated how accurately one can determine the binary parameters including wbd 
using advanced LIGO. Scharre and Will [23] and Will and Yunes [24] did the same analysis using LISA. Will and 
Yuncs improved the previous works by showing the constraint dependences on LISA position noise, acceleration noise 
and arm lengths. They also used slightly improved noise curve than the one used by Scharre and Will. However, these 
calculations do not include binary spins and also they used the pattern-averaged waveforms (having been averaged 
over the source direction {6s,(j)s) and the direction of the orbital angular momentum (^Lj^^l))- Berti et al. [25] 
calculated the determination errors including the effect of the spin-orbit coupling using LISA. Also, they performed 
the Monte Carlo simulations in which they randomly distribute the directions and the orientations of 10'* sources 
over the sky, calculate the constraint on wbd from each binary and take the average at the end. According to their 
analysis, for a (1.4-HlOOO)M0 a neutron star (NS)/black hole (BH) binary of SNR=\/200, 1 yr observation by LISA 
can put the bound wbd > 10799 on average. 

In massive graviton theories, the propagation speed of a gravitational wave depends on its frequency, which modifies 
the time of arrival from general relativity. This also affects the phasing of the gravitational waveforms. Will [26] 
included Xg into the binary parameters and carried out the matched filtering analysis, investigating how accurately 
one can determine Xg using advanced LIGO or LISA. Will and Yunes [24] did the same analysis using the improved 
noise curve for LISA. Again, they did not include the efl:cct of binary spins and also they used the pattern-averaged 
waveforms. Berti et al. [25] estimated the constraint on Xg by performing the Monte Carlo simulations including the 
additional parameters of the spin-orbit coupling (3 and the angles 9s,4's,(^l, and <j)L- According to their analysis, for a 
(10^ -I- 1O^)M0 BH/BH binary at 3Gpc, 1 yr observation by LISA can put the bound Xg > 1.33 x lO^^cm on average. 

In this paper, we improve the analysis of Ref. [25] by taking the following effects into account: (i) the spin-spin 
coupling, (ii) the eccentricity, (iii) the spin precession, (i) The spin-spin coupling a appears at the second post- 
Newtonian (2PN) order in the PN waveforms. Berti et al. reported that when they included both a and wbd or cr and 
Xg into parameters, they could not take the inverse of the Fisher matrix properly, because the ratio of the maximum 
eigenvalue to the minimum one of this matrix approaches the bound for the double precision computation. Hence, we 
carried out our analysis in the quadruple precision, (ii) As the binary radiates gravitational waves, it loses its energy 
and angular momentum [27]. Since this circularises the orbit as e oc f~^^^^^, usually the eccentricity is neglected in 
the analysis. We estimate how much the error of each parameter increases when we add eccentricity of 0.01 at 1 yr 
before coalescence into the binary parameters. The reason for choosing this (seemingly rather small) fiducial value for 
eccentricity is explained in Sec VII A. (iii) When the spin of the binary object is not zero, usually the spin-orbit eS'ect 
and the spin-spin effect cause the spin vectors Si (i=l,2) and the orbital angular momentum vector L to precess. 
The parameter estimation errors including effects of precession in general relativity have been calculated by several 
authors [28-33]. Vecchio [28] performed the Monte Carlo simulations of equal mass BH/BH binaries with LISA, taking 
only the leading spin-orbit terms in the precession equations. In this case, the precession equations can be solved 
analytically up to some appropriate order. This is called the simple precession approximation [34] . It has the property 
that L and Si precess around a fixed vector Jq which is almost the same as the total angular momentum vector J. 
Another property is that the following three quantities, LS, Si- S2 and S, become constant, where 5' = S'l + 5'2 is the 
total spin angular momentum and S is its magnitude. Lang and Hughes [29] solved the full precession equations, which 
include both spins of the binary stars, numerically and performed the Monte Carlo simulations of super massive black 
hole (SMBH)/SMBH binaries with LISA. For (10^ -t- 1O'^)M0 binaries at 2; = 1, they found that the determination 
accuracy for An/ij, increases by more than 2 orders of magnitude. Van der Sluys et al. [30-32] performed the Markov 
Chain Monte Carlo simulations for ground-based detectors, taking the simple precession approximation into account 
and performed the analysis that is beyond the Fisher matrix method. In Ref. [33] , they used the modeled gravitational 
wave signal injected into LIGO data so that the detector noise they used is more realistic than the Gaussian one. The 
precession gives some additional information to the waveform, which solves degeneracies among the parameters and 
enhances the parameter determination accuracy. Here, we estimate how much the effect of precession enhances the 
determination accuracy in alternative theories of gravity. We assume that the spin of one of the binary constituents 
can be neglected so that the analytic expressions in the simple precession approximation [28] can be applied. The 
so-called Thomas precession is neglected in Ref. [28] for simplicity but we found that this cannot be neglected for 
L-N » ±1, where N is the unit vector in the direction of the binary from the Sun. For this reason, we approximately 



3 



TABLE I: Comparison of the constraints on wbd and Ag in this paper and the ones obtained by several authors using LISA. 
Scharre and Will [23] , Will and Yunes [24] , and Will [26] performed the pattern-averaged analysis without including binary spins 
into parameters. Will and Yunes improved the previous works by showing the constraint dependences on LISA position noise, 

acceleration noise and arm lengths. They also used slightly better noise curve than the one used by Scharre and Will. Berti et 
al. [25] included spin-orbit coupling and performed Monte Carlo simulations. Wc extended their analysis by taking spin-spin 
coupling, spin precession and eccentricity into account. For the constraint on a>BD, wc show the one without eccentricity since 
we expect that the prior distribution A7e > reduces the results to the ones without including it. 
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take this Thomas precession into account. Neglecting the spin of a NS is justified from observations [35]. 

Following the analysis of Berti et al. [25], we assume that we detect gravitational waves from inspiral compact 
binaries from 1 yr before the coalescences with LISA. We carry out the following Monte Carlo simulations. We 
distribute 10^ binaries over the sky and estimate the determination accuracy of binary parameters including wbd and 
Xg for each binary using the Fisher matrix method. We take the average at the end. For the LISA noise curve, we 
use the same analytical approximation presented in Ref. [36]. 

We first performed our analysis using the pattern-averaged waveform. For the error estimation in Brans-Dicke 
theory, the inclusions of both cr and eccentricity into parameters reduce the determination accuracy by an order of 
magnitude. In particular, including eccentricity affects the estimation more than just including a. However, we found 
that if we impose the prior information on eccentricity, the constraint on wbd remains almost as stringent as the one 
without including eccentricity into parameters. For the analysis of massive graviton theories, the inclusion of these 
parameters only changes the results by a factor of a few. In this case, the inclusion of a affects more than the inclusion 
of eccentricity. Also, the prior distribution on eccentricity does not affect the result for the massive graviton theories. 
Next, we performed the Monte Carlo simulations including precession. We found that in the Brans-Dicke case, the 
results are not so much affected by taking precession into account. Using a NS/BH binary of (1.4 -I- 1OOO)M0 with 
SNR=v'200, estimation with all cr, eccentricity and precession taken into account leads to a constraint wbd > 2838 on 
average. When we include the prior information on eccentricity, we expect that the constraint becomes wrd > 6944, 
which is the same as the one without including eccentricity. In the case of massive graviton theories, the inclusion 
of precession has a more remarkable effect. Using a BH/BH binary of (10^ + 1Q^)Mq at 3Gpc, estimation with all 
cr, eccentricity and precession taken into account can constrain Xg > 3.06 x lO^^cm on average. This constraint is 
2.3 times stronger than the result obtained in Ref. [25]. To compare our results with the ones obtained earlier using 
LISA [23-26], we list them in Table 1. Since all of these analyses do not take into account the errors coming from 
the use of approximate waveforms and the limitation of Fisher analysis [38], the results shown in Table I are only 
approximate estimates. 

The organisation of this paper is as follows. In Sec. II, we review the waveforms of inspiral compact binaries in 

alternative theories of gravity. We discuss how df /dt is affected by the gravitational dipole radiation and the change 
in the propagation speed of gravitational waves. In Sec. Ill, we explain the output waveform of the detectors. They 
are the superposition of the two polarised waves. We consider the case of using two detectors and show the Fourier 
component of the restricted 2PN waveforms. Section IV discusses the effect of precession. In Sec. IV A, we review how 
the angular momentum vectors L and Si precess under the simple pre(;cssion approximation. In Sec. IV B, we discuss 
how the detector output expressed in Sec. Ill changes when the orbital angular momentum L evolves with time. We 
also show how we treat the so-called Thomas precession. In Sec. V, we describe how to estimate the determination 
errors of binary parameters using the matched filtering analysis. In Sec. VI, we present the noise curve of LISA that is 
needed when calculating the Fisher matrix. In Sec. VII, we explain the setups and present the results of our numerical 
calculations. In Sec. VII A 1, we use the pattern-averaged waveforms so that the angles 9s,4'S,0h,4'L are not taken 
into parameters. In Sec. VII A 2, we perform Monte Carlo simulations of 10^ binaries distributed over the sky. We 
calculate the determination errors for each binary and take the average. In Sec. VII B, we carry out the same Monte 
Carlo simulations including the effect of the simple precession approximation. In Sec. VIII, we summarise our main 
results and discuss some future works. 
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II. BINARY GRAVITATIONAL WAVEFORMS IN ALTERNATIVE THEORIES OF GRAVITY 



First, we calculate the gravitational waves coming from a binary system composed of compact objects in Brans- 
Dickc theory and massive graviton theories [25]. As the binary radiates gravitational waves, the orbit gets eirciilarised. 
Since we assume the observation starts 1 yr before coalescence, the eccentricity of this binary orbit is expected to be 
considerably small. Therefore, we only include the leading term for the eccentricity in the phasing part. 

In deriving the waveform, we use the post-Newtonian formalism, an expansion in terms of gravitational potential 
U and , where v is the typical source velocity. The emitted gravitational waves are the superposition of harmonics 
at multiples of the orbital frequency. The waveform h{t) can be written schematically as 



/i(t) = Re^ft;^(i)e'™^->=W , (1) 

\n,m / 

where n labels PN order and m is an index of the harmonics [40]. (/)orb — /* VL{t')dt' is the orbital phase, where ri(t) is 
the orbital angular velocity. In this paper we use the "restricted 2PN waveform." For the amplitude, we only take the 
leading Newtonian quadrupole term h^, and for the phase part, we use (t>orh{t) valid up to 2PN order. This is because 
the correlation between two waveforms is much more sensitive to the phase information than to the amplitude when 
we perform the matched filtering analysis. In this approximation the waveforms of + and x polarisations are 



h+{t) = A+ con (p{t), (2) 
h^{t) = A^sm(f>{t), (3) 



where (j){t) = 2^orb(i) and 



A. = '-^{i+it-m (4) 

A. = -'^(L-N). (5) 

Here mi, TO2 are the two masses, r is their orbital separation, Dl is the luminosity distance between the source and 
observer, L is the unit vector parallel to the orbital angular momentum, and TV is the unit vector pointing toward 
the source from the observer. Please see Appendix A for more explanations on polarisations. 

Next, we need to calculate the phase ^{t). First, we introduce the following useful mass parameters: 



M = mi + m2, (6) 

" - W («) 

M s ic'/^M"''' = M,r"^. (9) 

They represent the total mass, the reduced mass, the symmetric mass ratio, and the chirp mass, respectively. The 
rate at which the frequency changes because of the back-reaction due to the emission of gravitational waves is given 
by [25, 35, 41] 



^ ^ ^7r^/3_/V^5/3 jll/3 

dt 5 



, 157^ _i9,e 5 -1 96^ 3/5 /^743 
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24-^- ' 48- ' 5'^^'' l 336 ^T'^''' 



-I- (47r - /3)a;3/2 ^ 



/ 34103 13661 
V 18144 2016 



59 
18' 



r? + —nf +cr \ 



(10) 



where we defined the squared typical velocity, x = v'^ = {ttM J)^/^ . The first term in the square brackets represents 
the lowest order quadrupole approximation of general relativity. The second term is the contribution from small 
eccentricity [42]. 1^ = {ttMY'^/^ c^Jq^^'^ is the dimensionless asymptotic eccentricity invariant. (Note that our def- 
inition of Je is different from that in Ref. [42] by a factor (ttM)^^/^ so that we set our Jg dimensionless.) /o is an 
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arbitrary reference frequency and eg is the eccentricity when the frequency of quadrupole gravitational waves is /q. le 
is constant up to the leading order in e. The third term represents the dipole gravitational radiation in Brans-Dicke 
theory [22, 23]. u) = is the inverse of the Brans-Dicke parameter. <S = S2 — si with 



which is called the sensitivity of the i-th body. Here, GcS is the effective gravitational constant near the body and 
is proportional to the inverse of the Brans-Dicke scalar field there. The subscript denotes that we evaluate Sj at a 
sufficiently large separation from the body. This sensitivity roughly equals to the binding energy of the body per unit 
mass. For example, the sensitivities of white dwarfs and neutron stars are swd ~ 10~^ and sns ~ 0.2, respectively. 
Because of the no hair theorem, black holes cannot have scalar charges and hence sbh = 0.5. From Eq. (10), the 
contribution of dipole radiation becomes large as <S increases. Binaries with large S are the ones composed of bodies 
of different types. 

The fourth term is the contribution from the mass of the graviton [26]. When the graviton is massive, the propa- 
gation speed is slower than the speed of light, which modifies the gravitational wave phase from general relativity. Pg 
is given by [25, 26] 



xjii + zy ^'^^ 

where z is the cosmological redshift and Xg is the graviton Compton wavelength. For a flat Universe {^1^=0, ^a+^m = 
1), the distance D is defined as 



D.i±ir ^ (13) 

Ho Jo (l + z')Vf^M(l + 2')' + f^A 

Hq = 72km/s/Mpc is the Hubble constant, and flM — 0.3 and J^a — 0.7 represent density parameters of matter and 
dark energy, respectively. 

The remaining terms are the usual higher order PN terms in general relativity. The quantities j3 and a represent 
spin-orbit and spin-spin contributions to the phase, respectively, given by 



^ ^ uT.X^(^^^1± + 75v)L-S,, (14) 
i=i ^ ' 

a = ^xiX2(-247Si-S2 + 721(L-Si)(L-S2)), 

(15) 

where Si{i = 1,2) are unit vectors in the direction of the spin angular momenta. The spin angular momenta are given 

by Si = XiTT^i'^i where Xi are the dimensionless spin parameters. For black holes, they must be smaller than unity, 
and for neutron stars, they are generally much smaller than unity. It follows that |/3| < 9.4 and |c7| < 2.5 [25]. 

From d(j)/dt = 27r/, we get the following time and phase evolution of the gravitational radiation by integrating 
Eq. (10) with time, 



^ 157^ ig/6 52 , 128^ 2/5 4 /743 11 , 

8, ,/2 / 3058673 5429 617 2 

--(4.-^).3/^+2(^^^ + ^,+ -^2 



rj — a ] X 



(16) 



0(/) = <^c 
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(ttM/) 



-5/3 



785 , -.a/a 25 o 1 ^ 2/"^ 5 /743 11 
344^^^^ ^ - 336'^ - ^^^^^^^ + 3 (336 + ) - 

5/. «^ 3/2 ./ 3058673 5429 617 2 \ 2" 



(17) 



where tc and are the time and phase at coalescence. 
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III. DETECTOR RESPONSE 




FIG. 1: We use two types of coordinates: (i)a barred barycentric frame {x,y,z) tied to the ecliptic and centred in the solar 
system barycentre, (ii) an unbarred detector frame {x,y,z), centred in the barycentre of the triangle and attached to the 
detector. 



LISA is an all-sky monitor with a quadrupolar antenna pattern and consists of three drag-free spacecrafts arranged 
in an equilateral triangle, 5x lO^km apart. Each spacecraft contains a free- falling mirror so that LISA forms three-arm 
interferometers with opening angles 60°. The barycentre of each triangle orbits the Sun 20° behind the Earth and 
the detector plane is tilted by 60° with respect to the ecliptic. 

Following Ref. [43], we introduce two Cartesian reference frames: (i) a barred barycentric frame {x,y,z) tied to 
the ecliptic and centred in the solar system barycentre, with z (unit vector in z direction) normal to the ecliptic and 
xy-plane aligned with the ecliptic, (ii) an unbarred detector frame {x,y,z), centred in the barycentre of the triangle 
and attached to the detector, with z normal to the detector plane (see Fig. 1). The orbit of the detector barycentre 
can be written as 



9{t) = 7r/2, (j>{t) = 2TTt/T, (18) 

where T = lyear and we have assumed <^(< = 0) = 0. 

The interferometer having three arms corresponds to having two individual detectors. We label each detector as 
detector I and II, respectively. The waveforms measured by each detector are given as 



ha{t) = -^ ^^"^^ Apol,a(t) cos [(t){t) -f- tf-po\,a{t) + ^D{t)] , (19) 

where a = I, II represents the detector number (see Appendix A for more details). The polarisation amplitude 
^poi,a(Oi the polarisation phases </3poi,a(0 and the Doppler phase (fioit) are defined as 



ApoiMt) = \J{l + {L-NfYF+{tf+A{L-NfF^{t)\ (20) 

sin(^po,,„(t)) = ^^^.^'^f.^'K (22) 
'foit) - 27r/(t)i?sin^sCOs[<^(i)-^s], (23) 
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where L is the unit vector parallel to the orbital angular momentum and TV is the unit vector pointing toward the 
centre of mass of the binary, and and are the beam-pattern functions of + and x polarisation modes for each 
detector, shown in Appendix A. (^s- 0s) represents the direction of the source in the solar barycentre frame and R 
represents 1 AU. The Doppler phase denotes the difference between the phase of the wave front at the detector and 
the phase of the wavefront at the solar system barycentre. It arises from the motion of the detector around the Sun. 
The factor •\/3/2 in Eq. (19) comes from the 60° opening angle of adjacent detector arms of LISA. 

Later, we estimate the accuracy of determination of the binary parameters using the matched filtering analysis, 
where we work on the Fourier domain. Therefore we calculate the Fourier transform of the signal, 



(24) 



To evaluate this quantity, we use the stationary phase approximation [40]. When conditions d\nA/dt <C d(j)/dt and 
d^(l)/dt^ <C {dcp/dt)-^ are satisfied, the saddle point method can be used and the Fourier transform of a function 
B{t) = A{t) cos^{t) becomes 



B(/).l.,.,(| 



-1/2 



expi(27r/t(/)-0(/)-7r/4). 
Under this approximation, the Fourier component of the waveform h{f) becomes [25] 



(25) 



where the amplitude A and the phase ^'(/) are given by 



(26) 



A = 
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Also, when we integrate out the angle dependence from the waveform (26), it becomes 



(28) 



M/) = ^^r'/v*(/). 



(29) 



IV. PRECESSION 

In this section, we introduce an additional effect, the precession. The spin-orbit interaction and spin-spin interac- 
tion change the orientations of the orbital angular momentum vector L and the spin vectors Si. These vectors precess 
over a time scale longer than the orbital period but shorter than the observation time scale. This effect drastically 
changes the detected waveforms. 



A. Simple Precession 



In this paper, we assume that one of the spins of the binary constituents is negligible (i.e. Si 0). Then, the 
precession equations become 
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32 m' fM-^'''^ 



5 = 0, (31) 
i=f2+^^Uxi, (32) 



2 mi J 

(33) 



2 mi J 

where J is the total angular momcntiim J = L + S . Under this simple precession approximation, the following 
quantities become constant during the inspiral, Si ■ S2, k = L ■ S, and the magnitude of the total spin angular 
momentum S = |Si + S2I. Prom the above equations, it can be seen that both L and S precess around J with 
angular velocity 



In general, the precessing time scale is shorter than the inspiral time scale L/\L\. Therefore J changes in 

magnitude but the direction is almost constant. Then, the analytic form of L can be obtained up to some approximate 
orders (see Appendix B). 



B. Detector Response 

When precession is taken into account, the principal axis N x L varies with time so that the GW phase ^{t) no 
longer equals to twice the orbital phase ^orb(t) = / fl{t)dt. We define this difference 25(j){t) by 



(/.(t) = 2<^orb(t) + 2(50(0 . (35) 

This 54i{t) is the so-called Thomas precession phase. Wc specify the constant of integration so that (pothitc) = 4'c- In 
general, S(j){tc) 7^ so that (/>c does not equal to 0(ic) anymore in this case. From Eqs. (A4) and (35), the detector 
output ha{t) becomes 



hc.it) = ^A+F+{9s,h,i^s)cos2{^orh{t)+S<p{t)) + ^A^i^^^ (^s, ^s, V's) sin2(0orb(i) + 

= ^^^^(i^rcos20o,b + i^rsin2</.orb), (36) 
where a = I, II labels the detector number and F':"^ and F^™ are defined as follows: 



F™^ = {1 + {L ■ Nf)F+ cos2d(P - 2{L ■ N)F^ sm25(j), (37) 
Fi"" = -{l + {L-Nf)F+sm25(l>-2{L-N)F^cos26(t). (38) 

Following Eq. (19), we express this output (36) in terms of an amplitude and phase form, and we also take the motions 
of the detectors into account to obtain 



2 rD 

where A^'^^'^^it) and Vpo'i!a(*) ^^'^ gi^®'^ ^ 



hcit) = ^^^^X W cos(20orb(t) + CotJi) + (39) 



9 



A 




FIG. 2: The Thomas precession phase S<l){t) is the angle from the vector u to the principal axis N x L. We define 50' as the 
angle from the vector L x u to the one N x L. Notice that these two vectors always lie on the orbital plane. Sff)' equals to 
Sfj) + Tr/2 up to 0{\5L\/L). 



KolJt) = ViFT^W+iFFW, (40) 
cos«r«(t)) = (41) 

M^l'Zit)) = -^^y (42) 

The quantities {L ■ N),{L ■ ^) a-nd [N ■ [L x £)], which are needed to compute the polarisation angle 4's{i) in the 
beam-pattern coefficients F+ and , are expressed in Appendix B. 

Finally, we need to calculate the Thomas precession phase 5(t){t) . Apostolatos et al. [34] derived an explicit form as 



This expression includes an integration which makes the computational time very long. Vecchio [28] estimated the 
binary parameter accuracies by choosing a few random sources with and without including 5(j){t) and found that this 
term did not affect the results, and concluded that this term could be neglected. This is true for the binaries for which 
L ■ N never becomes close to ±1. However, when L ■ N ^ ±1, the direction of the principal axis (TV x L) changes 
rapidly with time, and hence the polarisation angle "03(0 also changes rapidly. Since it is the Thomas precession 
phase 5(j)(t) that cancels this rapid change, 5(j){t) cannot be neglected in this case. 

The direction which 5(t) is measured from is arbitrary. It can be seen from Eq. (43) that Apostolatos et al. [34] 
defined 5(j> as the angle measured from the principal axis at the time of coalescence. In general, since (TV x L)t=t^ 
does not lie on the orbital plane at a time t, we need to follow the evolution of the principal axis to calculate the 
Thomas precession S(j>{t). This is the reason why the integration appears in Eq. (43). 

Here, we try to derive an approximate expression of S(t)(t) that is not in the integral form. To do so, we use a 
specific vector L x u that always lies on the orbital plane, where m is a constant unit vector. In the limit of a large 
separation (i.e. t — >■ — oo or L — >• oo), L approaches to Jq. Using this Jo, we take u = N x Jq (see Fig. 2). We 
define an approximate Thomas precession phase 6(j){t) by the angle from this vector ii to the principal axis N x L. 
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(Wc may take the vector ii to be the principal axis at the time of coalescence (TV x L)t=t^- Then, S(f>(tc) becomes 
and (j){tc) = (t^c- However, the computation is easier for choosing w to be AT x Jq.) We denote the angle from 
-L X ti to iV X £- as 5(1)' , and the difference between Jq and L as 5L; L = Jq + SL. Then, 6^' can be related to S(j) 
as S(p' = 5(f) + 7r/2 up to 0{\5L\/ L). As 5(f>' always lie on the orbital plane, by using this angle we can express 5(l>{t) 
without integration, sin 50 and cos (50 can be written as 



sin 5(f) Pi — cos 5^' 

{LxN)-{Lx u) 
\L X 7V||i X u\ 
{L-N){L-u) 

^(l-(L.iV)2)(l-(L.n)2)' 
cos 6(f) « sin 50' 

N-{Lxu) 



Therefore the Thomas precession phase 5(f){t) can be expressed as 



(44) 



(45) 



50 Ri arctan 



{L ■ N){L ■ u) 



N-{Lxu) J 

The explicit forms of L ■ N, L ■ u and N ■ {L x u) are given in Appendix B. 



(46) 



V. PARAMETER ESTIMATION 

We use the matched filtering analysis to estimate the determination errors of the binary parameters 6 [40, 44]. 
Wc assume that the detector noise is stationary and Gaussian. "Stationary" means that different Fourier components 
n{f) of the noise are uncorrelated and we have 



{n*{f)n{n)=S{f-fhSM)- 



(47) 



Here, (...) denotes the expectation value and Sn{f) is the noise spectral density. Then the noise takes the Gaussian 
probability distribution given by 



p{no) (X exp 



df 



\MfW 



(l/2)5„(/) 



oc exp 



(48) 



where we have defined the inner product as 



{A\B) = Re 



= 4Re 



df 



•oo 



Jo 



A*{f)B{f) 



(49) 



(/) 



The signal to noise ratio(SNR) for a given h is 



p[h] = Vihlhj. 



(50) 
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The detected signal s{t) = h{t, 9i) + no(t) is the sum of the gravitational wave signal h{t; Ot) and the noise n{t), 
where Ot is the true binary parameters and no(t) is the noise of this specific measurement. Then, Eq. (48) becomes 



{ht\s)-\{ht\ht) 



(51) 



where p^^\Ot) represents the distribution of prior information and ht = h{6t). We determine the binary parameters 
as 6 that maximise the probability distribution p{6t\s). Then d is the solution of the following equation, 



{diht\s)-idiht\ht) = o, 



(52) 



where di = We denote the error in the determination of 0^ as A^'; 9^ = 6^ + AO^. Next, we expand Eq. (51) in 
powers of A^' up to quadratic order and we get 



p(6>|s)ap(°)(6')exp 



(53) 



where Tij = {didjh\h — s) + {dih\djh) is called Fisher matrix. Since h — s = —n, in the limit of large SNR, we can 
neglect the first term of Tij and 



Tij = {dih\djh). 

When we use both detectors I and II, the Fisher matrix becomes 



(54) 



Tij = {dihi\djhi) + {dihu\djhu). 
We take into account our prior information on the maximum spin by assuming 



(55) 



p^°\e) (X exp 



-(-Y --(-)' 

2 V 9.47 2 V2.5/ 



(56) 



Then, the rms of A9^ can be calculated by taking the square root of the diagonal elements of the covariance matrix 
E*-' , which is the inverse of the Fisher matrix: 



and Fij is defined by 



S*^' = {AO' AO') = [f-^f 



(57) 



p(°)(6>)exp 



--ViiAe'Ae^ 



= exp 



-hijAe'AO^ 



(58) 



VI. LISA NOISE SPECTRUM 



In this section, we introduce the noise spectrum of LISA following Ref. [36]. The instrumental noise spectral 
density for LISA is given as 



9.18 X 10"52 



/ 
1 Hz 



-4 



1.59 X 10-41 + 9.18 X 10-38 ( -1- 

\ IHz 



Hz" 



(59) 
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FIG. 3: The noise spectral density for LISA. 



Besides instrumental noise, there are foreground confusion noises. These confusion noise spectral densities and the 
energy densities of gravitational waves are related as 



where pc = is the critical energy density of the Universe and 



(60) 



r^Gw = 



1 rfpGW 

Pc dlnf 



(61) 



is the energy density of gravitational waves per log frequency normalised by pc- The energy density of gravitational 
waves coming from extra-galactic white dwarf binaries has been estimated as 51gw — 3.6 x 10~i^'^(//10^^Hz)^/'^ and 
the noise spectral density becomes [45] 



ST~^"\f) = 4.2 X 10-4^ 



1 Hz 



-7/3 



Hz-1. 



(62) 



On the other hand, the energy density of gravitational waves coming from galactic white dwarf binaries has been 
calculated as 50 times larger than the one coming from extra-galactic white dwarf binaries. Therefore the noise 
spectral density becomes [46] 



Sfif) = 2.1 X 10-45 



1 Hz 



-7/3 



Hz-1. 



(63) 



We compute the total noise spectral density as 



Shif) = min 



qmst f -f\ 



^cx— gal 



(/)• 



(64) 



Here dN/ df is the number density of galactic white dwarf binaries per unit frequency, for which we use the estimate 
given in Ref. [47] 



^ . 2 X 10-3HZ- 
df V 1 Hz 



-11/3 



(65) 



K ~ 4.5 is the average number of frequency bins that are lost when each galactic binary is fitted out. This noise curve 
is drawn in Fig. 3 
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VII. NUMERICAL CALCULATIONS AND RESULTS 

Following Berti et al. [25] , we estimate how accurately we can determine binary parameters by detecting inspiral 
gravitational waves with LISA. The total noise spectral density for LISA is given by Eq. (64) and we introduce cutoff 
frequencies for LISA as follows: 

/low = lO-^Hz, /high = IHz. (66) 

Here, wo use the matched filtering analysis, assuming that the observation starts 1 yr before the coalescence. We 
assume that we use two detectors which corresponds to one triangle interferometer as mentioned in Sec. III. We 
numerically calculate the determination error of each parameter which is given by the square root of the corresponding 
diagonal element of S*-' . We use the Gauss- Jordan elimination [48] for inverting the Fisher matrix Fy to obtain S*-' . 
In calculating this r^ , we need to perform the integration in the frequency domain. Following Ref. [25], we take the 
frequency range of this integration as (/in, /fin); 

/in = max{/iow, Ayr}, /fin = min{/high, /isco}, (67) 

where 

/lyr = 4.149 X 10"^ 
is the frequency 1 yr before the coalescence, and /isco is given by 

= 63/27rMBH- ^^^^ 

For the waveforms, we use the restricted 2PN waveforms. We only take up to 2PN because spin terms are known 
only to this order. There are 15 parameters in total: the chirp mass InA^, the dimensionless mass parameter Inry, 
the coalescence time tc^ the coalescence phase 0c, the distance to the source InD, the spin-orbit coupling coefficient 
13. the spin-spin coupling coefficient cr, the inner product of the orbital angular momentum and the total spin angular 
momentum k, the precession angle ac that characterises a at the time of coalescence (the precession angle a is 
defined in Appendix B), the dimensionless asymptotic eccentricity invariant /g, (^S, <?^s) for the initial direction of the 
source, (0,j, (/)j) for the initial direction of the total angular momentum, and finally, w, the inverse of the Brans-Dicke 
parameter, or j3g, the quantity that is proportional to the square inverse of the graviton Compton wavelength A^, 
depending on which theory we are aiming to constrain. For constraining the Brans-Dicke parameter, we consider 
NS/intermediate mass black hole (IMBH) binaries, assuming the difference between NS and BH sensitivities S to be 
0.3. For constraining the graviton mass, we consider SMBH/BH binaries at 3Gpc. In this case, we fix the cosmological 
parameters as follows: fi^ = firad = 0, = 0.3, = 0.7. We take the Hubble constant to be Hq = 72 km/s/Mpc. 
There are three main differences from Ref. [25]. 

(i) We include the spin-spin interaction a: Berti et al. [25] reported that when they included both a and the 
parameter for alternative theories of gravity like Cu, the ratio (we denote this by R) between the smallest and the largest 
eigen values of the Fisher matrix T approached their machine's floating-point precision, and they could not obtain 
enough accuracy for the inversion of this matrix. We evade this problem by performing our numerical calculations in 
quadruple precisions. Also, we use the trick explained in Appendix C to make sure the numerical inversion is correctly 
performed. Basically, we rcscalc the basis vectors so as to make all the diagonal components of the Fisher matrix 
equal to 1 [36] , then take the inverse of this normalised matrix, and finally multiply the factor for the rescaling back 
to obtain the inverse of our original Fisher matrix. Even if the ratio R for the original Fisher matrix approaches the 
machine's floating-point precision, the one for the normalised Fisher matrix stays smaller than the floating-point. To 
check that our matrix inversion has been performed correctly, we followed a similar procedure described in Berti et 
al. [25]. We also mention this in Appendix C in more detail. 

(ii) We take the eccentricity into account: The binary system gradually loses energy and angular momentum 
becaiisc! of the gravitational radiation which circularises the orbit. In fact, since the eccentricity decreases following 
e (X [27], usually one does not include the eccentricity into binary parameters assuming that it is a priori. 
Still, we should take the eccentricity into account for more practical analysis. 

(iii) We also take the spin precession into account: In this case, the orbital angular momentum L oscillates, 
which introduces the additional oscillations in both the amplitude and the phase of the waveform. This additional 



M 



IO6M0 



-5/8 



-^obs 

lyr 



-3/8 



(68) 
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information solves degeneracy in the binary parameters, enhancing the determination accuracy. The parameter 
estimation including precession has been estimated by several authors using LISA [28, 29] and using the detectors 
on the ground [30-33], in the framework of general relativity. Our calculation is the first one to include precession 
in parameter estimation in the context of modified gravity. We consider the simple precession approximation by 
assuming that one of the spins of the binary constituents is 0. References [28, 30-33] estimated the determination 
errors of the binary parameters under this simple precession approximation. 

A. No precession 

In this subsection, we show the results without spin precession effects. In this case, k and ac are excluded from 
our binary parameters. Following Ref. [25], we perform the numerical integration using the Gauss-Lcgcndrc routine 
GAULEG [48]. Gauss-Legendre quadrature takes the abscissas at the zeros of the A^-th Legendre polynomials. With 
this method, integrand polynomials up to (2^" — l)-th order can be calculated exactly. We take the number of 
frequency bins to be 600. The fiducial values of the parameters are tc = ipc = P = cr = u) = fig = Q. 

In Ref. [49], Hopman and Alexander performed Monte Carlo simulations in which they followed a star on a relativistic 
orbit and added small perturbations to simulate energy dissipation and random two-body scattering. For (1.4-|-1O'^)M0 
NS/BH binaries, the probability distribution of eccentricity at / = 2 x 10~''Hz peaks at e = 0.998. From this 
condition, we further evolve the orbit to yield e = 0.026 at 1 yr before coalescence. (See Appendix D for more details 
on calculations of the eccentricity evolution.) For SMBH binaries, Armitage and Natarajan [50] have simulated the 
interaction between a binary and its surrounding circumbinary gas disk. For M = lO^M© SMBH binaries, the 
eccentricity at 1 yr before coalescence is e = 0.02 for the equal mass binaries and e = 0.04 for the unequal mass 
binaries with mass ratio 0.1. Recently, Bcrcntzcn et al. [51] performed A'^-body simulations of SMBH binaries in 
rotating galactic nuclei. They followed the evolution from kiloparsec separations to the final relativistic coalescence, 
including post-Newtonian corrections up to 2.5. They found that when SMBH binary reaches the separation of 100 
Schwarzschild radii, the typical eccentricity becomes e = 0.04. When the total binary mass is M = lO^Mo, the 
gravitational wave frequency at r = 100 Schwarzschild radii is / = \J M / {ir'^r^) = 2.25 x 10~^Hz. This frequency is 
roughly the same as the one at 1 yr before coalescence. Taking these results into account, we set the fiducial values 
of eo as eo = 0.01 at /o = /ly,- for both NS/BH and SMBH binaries. 

For the estimations without precession, we assume that the fiducial values of the binary spins arc zero. In these cases, 
the total angular momentum J is equivalent to the orbital angular momentum L. Therefore we take the direction 
of the orbital angular momentum (^l,</'l) as binary parameters instead of {Oj,(j)j). We treat the angles {0s,(l>s) and 
iOL,(pL) in two different ways. First, we take the average of these angles and calculate the determination accuracy 
of the binary parameters. We call this pattern-averaged estimate and we use only 1 detcc;tor for the analysis. The 
second one is the Monte Carlo simulation. We randomly distribute 10^ binaries (10^ sets of these angles), calculate 
the determination accuracy for each binary, and take the average at the end. In this case, we use 2 detectors for the 
analysis. 

1. Pattern-averaged estimates 

In the case of pattern-averaged estimates, there are only 9 parameters in total: In A^, Inrj, tc, ^c, Dl, /3, cr, 1^ and 
u) or pg. The waveform is given by Eq. (29). The derivative of this waveform with respect to each parameter is taken 

analytically as shown in Appendix E. 

Brans-Dicke theory.- 

Here, we consider the parameter estimation in the context of Brans-Dicke theory. We think of four NS/BH 

binaries of SNR=10. When performing Fisher matrix analysis, this SNR might be too small [52] and we should use 
the binaries with much larger SNR. In that case, since the constraint on wbd is proportional to SNR, all we have to 
do is to appropriately scale the results obtained for SNR=10. We fix tons = 1.4M0 and take rriBH = 400, 1000, 5000 
and 10^ M0. For each binary, we show three results in Table II: the first line represents the determination accuracies 
of binary parameters without taking a and /e into parameters (the same estimates by Berti et al. [25]), the second 
line shows the results including a but not 1^ into parameters, and the third line shows the ones including both a and 
/e into parameters. 

Prom this table, one can see that adding parameters reduces the determination accuracy. This is because the 
parameters are strongly correlated and adding parameters dilutes the binary information in the detected gravitational 
waves. Including both a and into parameters increases the determination errors by roughly 1 order of magnitude. 
In particular, including 7e increases the errors on wbd more than just including a. The reason for this is as follows. 
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TABLE II: The results of error estimation in Brans-Dicke theory. These are calculated with pattern-averaged analysis for 
NS/BH binaries with BH masses 400Mo, IOOOMq, 5000Mq and IOOOOM0 for p = 10. We fix NS masses to IAMq. For each 
binary, the first line represents the results shown in [25] which does not include a nor 7e into parameters. The second line 

represents the results taking a into account, and the third line shows the ones that include both cr and /e- Here /g is the value 
of at / = 0.3Hz. We used only one detector for the analyses. 
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FIG. 4: The curve showing 68% confidence level on the Ae^-Aw plane for a (1.4 + 4OO)M0 NS/BH binary with p = 10, such 
that the fiducial values lie at the centre of the ellipse. 

In the phase \E'(/), the term containing wbd is proportional to f~^^^, whilst the ones containing a and Ig have the 
frequency dependence /*/^ and f~^^^^, respectively. The terms containing wbd and le both have negative power-law 
indices and the dependences on the frequency are similar. Therefore wbd has stronger correlation with le than cr. 

Comparing the constraints from the binaries having different BH mass, one can see that the constraint becomes 
more stringent as the BH mass decreases. This is because the orbital velocity of the binaries become slower, which 
makes the dipolc contribution relatively greater. 

In Fig. 4, we plot the curve of constant probability on the A/g-Aa) plane for a (1.4 -I- 4OO)M0 NS/BH binary with 
p — 10. The fiducial values lie at the centre of the presented ellipse at 68% confidence level. Since 1^ is always 
positive, it can be understood from this figure that we can decrease the upper bound on Aw by imposing the prior 
distribution le > 0. In fact, the upper limit for Au> is obtained at A/g = in this case. This means that when we 
include the above prior information, the constraint on lobu including eccentricity should be similar to the case without 
taking eccentricity into account. Therefore we conclude that the constraint on wbd is not so much affected by taking 
eccentricity into account. 

Massive Graviton theories.- 
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TABLE III: The results of error estimation in massive graviton theories. These are calculated with pattern-averaged analysis 
for SMBH/BH binaries with masses (lO'' + 10^)Mq, (10^ + 10®)Mq, (10® + W^)Mq and (10*^ + W^)Mq at 3Gpc. As in the 
Brans-Dicke case, the first line of each binary represents the estimation without taking a nor le into parameters. These results 
do not exactly match with the ones shown in [25] because we also considered prior information. The meaning of second and 
third lines are the same as in Table II. 
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Next, we eonsider the parameter estimation in the context of massive graviton theories. We also consider four 
BH/BH binaries in this case: (10^ + 10'^)Mq, (10^ + 10^)Mq, {10^ + 10^)Mq and (10^ + 10^)Mq. As in the Brans- 
Dicke case, we perform three types of analyses whose results are shown in Table III. Berti et al. [25] reported that the 
effect of adding prior information on the maximum spin is negligible. We confirmed that including prior information 
changes the results just about a few percent. Berti et al. have also claimed that when they include both a and w or (3g 
into parameters, they cannot take the inverse of the Fisher matrix properly. This is true but, if we include the prior 
information, the matrix inverse can be performed successfully. Hence, wc decided to include the prior information. 
This is the reason why the first column of each binaries in Table III slightly differs from the results shown in Ref. [25] . 




FIG. 5: The curve showing 68% confidence level on the Ae^-AjBg plane for a (10^ + 10^)Mo SMBH/BH binary at 3Gpc. The 
fiducial values lie at the centre of the ellipse as in Fig. 4. 

From Table III, it can be seen that in this case, including both a and /(. into parameters reduces the accuracy 
only by a factor of a few. This indicates that /3g is only weakly correlated with other parameters. In the massive 
graviton case, including a affects the constraint on \g more than including If.. The reason is the same as in the case 
of Brans-Dicke. Since both terms containing Xg and a in the phase \E'(/) have frequency dependences with positive 
power-law indices, their correlation is stronger than the one between Xg and le ■ 
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TABLE IV: The results of error estimation in Brans-Dicke theory for (1.4+1OOO)M0 NS/BH binaries without pattern averaging. 
We used two detectors for the analyses and we fixed p = V200 {p = 10 for each detector). We performed the following Monte 
Carlo simulations. We distribute 10* binaries, calculate the error of each parameter for each binary and take the average. The 
first half of the table shows the results without taking preeessional effect, and the second half represents the ones including 
precession. The first line of (^aeli part shows the ix^sults without taking I, into parameters, whilst in the second line this is 
taken into accouut. rr is included iu the paramotor tor all the cases. 
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Figure 5 shows the curve of constant probability on the A/e-A/3g plane for a (lO'^ + lO^)Af0 SMBH/BH binary 
at 3Gpc. The fiducial values lie at the centre of the ellipse at 68% confidence level. Compared to Fig. 4 of the 
Brans-Dicke theory, it can be seen that the ellipse is slightly tilted in the other way. This indicates that the constraint 
on Xg would not change even if we impose the prior distribution A/,. > 0. However, this ellipse is tilted only slightly, 
which suggests that /3g and Ig are correlated only weakly, as we have stated above. Therefore, the constraint on Xg 
is affected only weakly by the inclusion of into parameters. 

When the graviton has a finite mass, its phase velocity is expressed as [26] 



Therefore, at a lower frequency the difference between Vph and the speed of light (c = 1) is larger. This is why the 
heavier binaries put a more stringent constraint on Xg. 



2. Estimates without pattern-averaging 

In the case of estimates without pattern-averaging, there are 13 parameters in total: the 9 parameters appeared 
in the last subsection and 4 angles representing the source direction and the inclination of the orbit. We randomly 
generate 10'' sets of the four angles so that cos^s and cos6'l are uniformly distributed in the range [—1, 1] and (f)s and 
(jj-L in the range; [0,27r]. We use the RANI routine [48] for generating random numbers and calculate the parameter 
estimation errors for each binary. We take the average at the end. 

We take the parameter derivatives of the part e**^'^^ in the waveform analytically as in the pattern-averaged case, 
whilst we take the derivatives of the rest numerically. The angular resolution AOs is defined as 



AQs = 27r|sin^-s|ys;^;S^~^s|~. (71) 

We consider the following two cases: including not le but a into binary parameters, and including both a and 1^ into 
parameters. 

Following Ref. [25], we display the results in histograms. For each binary parameter, we obtain 10^ error values. 
We group them into A^bins bins as follows; if an error value X satisfies 



In(Xmin) + 



-^bins 



< \n{X) < 



In(Xmin) 



max 



(72) 



then we define that this belongs to the j-th bin. We divide the number of binaries of each bin by the total number of 
binaries 10000 to get a "probability distribution." We show this probability against the determination error of each 
parameter in a histogram below. 



Brans-Dicke theory.- 
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FIG. 6: The histograms showing the probabihty distribution of the lower bound of wbd obtained from the Monte Carlo 
simulations with (1.4+ 1000) Mq NS/BH binaries with p = \/200 in Brans-Dicke theory. The (black) dotted thin one represents 
the estimate without precession and 7e is not taken into account. The (blue) solid thin one shows the one without precession 
but 7e is taken into parameters. The (purple) dotted thick one includes precession but does not include 7e. The (red) solid 
thick curve shows the one including both precession and Je. Here the prior information Je > is not used to estimate the lower 
bound of ojbd 



We fix the masses of the binary constituents as mNS — IAMq and meH = IQ'^Mq. We consider binaries of 
p — VSOO (roughly speaking this corresponds to p = 10 for a single detector). We show the results in the upper half 
of Table IV. The first row represents the error estimation without including /g as a fitting parameter and the second 
row shows the one including /g. The corresponding histograms are shown in Figs. 6, 7 and 8. The (black) dotted 
thin line represents the results without including Jg as a binary parameter and the (blue) solid thin line shows the 
one including /g. The (purple) dotted thick and the (red) solid thick lines are the ones including precession, which 
are to be explained later. For the moment, let us focus on the (black) dotted thin and the (blue) solid thin lines only. 
Figure 6 is the histogram of the lower bound of wbd- Including /g increases the error by a factor of 4, as expected 
from the results of pattern-averaged estimation. Figure 7 shows the histograms of AAl/A^, A77/77, A/?, and Act. 
These accuracies also get worse by several factors when Ig are taken into account. Figure 8 shows the histograms 
of /S.Dl/ Dl and Ails. The (black) dotted thin and (blue) solid thin lines in the histogram of /S.Dl/ Dl are almost 
the same. As appears only in the amplitude of the waveform, this parameter is uncorrelated with the phase 
parameters such as /g. The histograms of /S.Dl/Dl have tails on the lower accuracy side. These tails are due to the 
binaries with L ■ N = ±1 because then the distance and the direction of L degenerate. The estimation errors for 
ADl/Dl go up to 10"' for our calculation but we only show them up to 10 in Fig. 8. Again, the accuracy of Afis 
goes down by several factors when wc include /g. 

In Figs. 6 and 7, there are also some tails in the histograms. The ones with worse accuracy are due to the binaries 
with L ■ N = 0. Figure 9 shows 10^ plots of the constraints on wbd against L ■ N, taking eccentricity into account 
(these plots correspond to the (blue) solid thin histogram in Fig. 6). It can be seen that the constraints get worse 
when the binaries become L ■ N = 0. The constraints on wbd against 6*3 and (j)s are shown in Fig. 10. It is clear that 
binaries with some special source directions {Os,(t>s) have the parameter estimation accuracy enhanced. These tails 
do not depend on the directions of L. We found that they are due to the motion of the detectors (the Doppler phase 
ipD{t)). The Eq. (23) for foit) can be recast into the form [53] 



f(t) _ _ _ 
^D{t) = 2tt^ sin 0s cos[</.(t) - 0s], (73) 

Jc 



with the critical frequency fc defined by 
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FIG. 7: The histograms showing the probability distribution of the error estimates of the chirp mass ^Mj M, the reduced 
mass the spin-orbit coupling A/3, and the spin-spin coupling Act in Brans-Dicke theory. These are obtained from the 

Monte Carlo simulations with (1.4 + 1OOO)M0 NS/BH binaries with p = \/200. The meaning of different types of curves is the 
same as in Fig. 6. 




FIG. 8: The histograms showing the probability distribution of the error estimates of the luminosity distance ADl/Dl and 
the angular resolution Af2g in Brans-Dicke theory. These are obtained from the Monte Carlo simulations with (1.4+ 1OOO)M0 
NS/BH binaries with p — \/200. The meaning of different types of curves is the same as in Fig. 6. 



/, = ^ - 2.00niHz, (74) 
it 

for R — lAU. This shows that the effect of ifioit) is important for frequencies higher than 2mHz. Since the relevant 
frequency range for a (1.4 -I- 1000)Mq binary is 3.66 x 10~^ — l.OOHz, (fioit) has some remarkable contribution. In 
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L-N 

FIG. 9: 10* plots of the lower bounds of wbd against L ■ N. This result is obtained from the Monte Carlo simulations with 
(1.4 + 1000)Mo NS/BH binaries with p = \/200 in Brans-Dicke theory. 7e is taken into account though the precessional effect 
is not included. 
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FIG. 10: 10* plots of the lower bounds of wbd against 9s and (/>s, as in Fig. 9. 



Figs. 11 and 12, we show the same plots but without including the Doppler phase. It is clear that the anomalous 
peaks have disappeared in this case. The comparison of the constraint on ojbd between these cases suggest that the 
degeneracies among some parameters become strong when we include (pD- However, these degeneracies are solved in 
some special cases which correspond to the peaks in Fig. 10. For example, the reason why = f is special can be 
understood as follows. Since (poit) is proportional to sin^s, the derivative of (fiD with respect to is proportional 
to cos^s- Therefore when = ^, this term vanishes and there would be no degeneracy between wbd and that is 
caused by the Doppler phase. 

Massive Graviton theories. - 

We fix the masses of the binary constituents to lO^M© and 10^ Mq, and the distance Dl =3Gpc. The results are 
shown in the upper half of Table V. As in the Brans-Dicke case, the first row represents the error estimation without 
including /g into binary parameters and the second row shows the one including /g. As in the Brans-Dicke case, 
the corresponding histograms are shown in Fig. 13 for the lower bound of Xg, in Fig. 14 for the error estimations of 
masses and spins, and in Fig. 15 for the ones of the distances and the source directions. Again, the (black) dotted 
thin line represents the results without including Ig as a binary parameter and the (blue) solid thin line shows the 
ones including I,,. In the massive graviton case, the parameter estimation is only weakly dependent on the inclusion 
of parameter Ig. 

There are 5 parameters in the phase ^{f) shown in Eq. (28); A4,r],P,a, and I3g. On the other hand, there are 

only 4 PN terms in total (the leading quadrupole term, IPN, 1.5PN, and 2PN). This means that these 5 parameters 
degenerate when we do not include precession. Then, one of the parameters is totally unconstrained without prior 
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FIG. 11: The same plot as Fig. 9 except we do not include the Doppler phase foit) in this case. 




distribution. The dispersion that the histograms for Xg, /U, ^ and a have in Figs. 13 and 14 is mainly determined 

by the level of uncertainty in the determination of the spin parameter a. When we do not take into account the 
precession, the level of uncertainty of a is totally determined by the prior distribution, and hence it is identical for all 
binaries. Therefore each histogram for Xg, fj,, (3 and a has a very sharp peak when we do not include the precession 
effect. 



B. Including Precession 



In this subsection, wo show the results when we include the spin precession. There are 15 parameters in total. 
We perform the Monte Carlo simulations as in the previous subsection. Here, we set the dimensionless spin parameter 
of the lighter body of binaries to and that of the heavier body to 0.5. We randomly distri bute the fiducial values of 
K = Z, • S' in the range [-1,1] and choose etc randomly in the range [0, 27t]. Wc use L — ii\J Mr for the calculation of 
orbital angular momentum, where the separation is given as r = M^/^/ (tt/)^/'^. (This is derived from = ^JM/r = 
(ttM/)^/-^.) When we include the spin precession, the waveform gets some additional oscillations. This worsens the 
precision of the polynomial approximation to the waveform. For this reason, we decide not to use the Gauss-Legendre 
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FIG. 13: The histograms showing the probabihty distribution of the lower bound of Ag obtained from the Monte Carlo 
simulations with (10^ + 10^)Mq BH/BH binaries at 3Gpc in massive graviton theories. The (black) dotted thin one represents 
the estimate without precession and le is not taken into account. The (blue) solid thick one shows the one without precession 
but Je is taken into parameters. The (purple) dotted thin one includes precession but does not take 7e into account. The (red) 
solid thick one shows the one including precession and le is also taken into account. 
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FIG. 14: The histograms showing the probabihty distribution of the error estimates of the chirp mass AA4/A4, the reduced 
mass A/x//x, the spin-orbit coupling A/? and the spin-spin coupling Act in massive graviton theories. These are obtained from 
the Monte Carlo simulations with (lO'^ + 10®)Mq BH/BH binaries at SGpc. The meaning of different types of curves is the 
same as in Fig. 13. 
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TABLE V: The results of the Monte Carlo simulations in massive graviton theories for (10^ + W^)Mq BH/BH binaries at 3Gpc 
without pattern averaging. We used two detectors for the analyses. As in the Brans-Dicke case, we distribute 10^ binaries, 
calculate the error of each parameter for each binary and take the average. The meaning of each row is the same as in Table IV. 



Cases 




SNR 


AlnM(%) 


Alnr; 


A/J 


AIu^Dl 


Ans(lO-^str) 


Aa 


No precession 


















Excluding 1^ 


0.40598 


1540 


0.0507 


0.841 


6.27 


0.0230 


0.957 


1.89 


Including le 


0.40507 


1540 


0.191 


0.927 


6.54 


0.0233 


0.972 


1.94 


Including precession 

Excluding le 
Including le 


4.8540 
3.0570 


1596 
1586 


0.00838 
0.0269 


0.00675 
0.00708 


0.0117 
0.0120 


0.00189 
0.00192 


0.366 
0.364 


0.0508 
0.0825 




FIG. 15: The histograms showing the probability distribution of the error estimates of the luminosity distance ADl/Dl and the 
angular resolution AQs in massive graviton theories. These are obtained from the Monte Carlo simulations with (10'^ + 10^)Mq 
BH/BH binaries at 3Gpc. The meaning of different types of curves is the same as in Fig. 13. The accuracy of the luminosity 
distance with and without precession are shown separately. 



quadrature for the numerical integrations. Instead, we cut the integrand into 10^ pieces, equally binned in terms of 

Brans-Dicke theory.- 

We fix the masses of the binary constituents to rrijsss — 1.4M0 and msH = lO^Mo and SNR p to \/200. We neglect 
the spins for the neutron stars, which is supported from observations [35]. For the dimensionless spin parameters of 
the black holes, we adopt x = 0.5. The results are shown in the lower half of Table IV. The first row in this table 
represents the results without taking eccentricity into parameters and the second row shows the results including 
Ig. The corresponding histograms are again shown in Figs. 6, 7 and 8. The (purple) dotted thick line represents the 
estimation without taking I^, into parameters and the (red) solid thick one is the estimation with Jg. Figure 6 shows 
that the constraints on wbd increase by 20% when the precession effect is taken into account. In [25], Berti et al. 
found that by detecting (1.4 + 1OOO)M0 NS/BH binary gravitational waves of p = -\/200 with LISA, one can put a 
constraint cjbd > 10799 on average. Comparing this result with ours, it is understood that including the effects of cr, 
le and precession, in total, leads to weaken the constraint by an order of magnitude. Therefore these effects cannot be 
neglected. However, as we stated before, imposing the prior distribution A/g > would reduce the upper limit of Aui 
to the case without including /g into parameters. This means that this prior information strengthens the constraint 
on ujbd by a factor of 6. In this case, the constraint on wbd becomes wbd > 6944, which is just 1.6 times lower 
than the one obtained in [25]. Figure 6 shows that when eccentricity is taken into account, the effect of precession is 
larger than in the case without including it. This is because the degeneracy between parameters are disentangled by 
precession. In Fig. 7, one can see that the lower accuracy tails disappear when we include precession, as there is no 
binary with L ■ N fixed to due to the precession of L. Also in Fig. 8, the tails on the histograms disappear when 
we include precession for the same reason. 

We have checked that when we take the limit k = 1, our Monte Carlo simulation reduces to the one without 
including precession. However, we found that when we set k = 0.999, the constraint on wbd differs from the one 
with K = 1 by several factors. This unphysical fact shows that we cannot estimate the errors properly with Fisher 
analysis in some regions. Above k = 0.999 (and below k = —0.999), the estimation error of precession angle Aa (not 



24 



0.07 



fO 



;: ^BH=0.5 



006 



0.05 



prec. X 




0.01 



I 




1 




0.2 



05 



A/3 



FIG. 16: The histograms showing the probability distribution of the error estimates of the spin-orbit coupling A/3 in Brans- 
Dicke theory. These are obtained from the Monte Carlo simulations with (1.4 + 1000) Mq NS/BH binaries with p — \/200. 
Eccentricity is not included into the binary parameters. The (black) dotted thin one represents the estimate without spin 
precession and the fiducial value for xbh is set to 0, whilst the (black) solid one represents the one with xbh = 0.5. The 
(purple) dotted thick one represents the one including spin precession and the fiducial value for xbh is set to 0.5. 

Aoic) exceeds 27r. This means that the effect of precession is too weak to be determined and the linear analysis of 
Fisher matrix is invalid. Therefore we have to say that our calculation cannot be trusted for these regions. When we 
exclude these regions and recalculate the average for the constraint on wbd, it becomes wbd > 3252 when we include 
eccentricity. Comparing this with the corresponding result in Table IV, wbd > 3251, we can say that the contribution 
of these untrustable regions are not so important. 

In Fig. 7, the results of A/3 may look a bit strange. When le is not included into parameters, the determination 
accuracy gets worse when we include spin precession. This is due to the difference of the fiducial values that we 
are taking. When we do not consider the effect of spin precession, the fiducial value of the dimensionless BH spin 
parameter xbh is set to 0. On the other hand, when we include spin precession, it is set to 0.5. If we set this 
BH spin parameter to 0.5 and perform the Monte Carlo simulation without including precession, the histogram gets 
shifted. This is shown in Fig. 16. The (black) dotted thin one represents the estimate without spin precession and 
the fiducial value for xbh is set to 0, whilst the (black) solid one represents the one with xbh = 0.5 and k = 1. The 
(purple) dotted thick one represents the one including spin precession and the fiducial value for xbh is set to 0.5. 
Comparing the (black) solid histogram and the (purple) dotted thick one, it can be seen that the inclusion of spin 
precession enhances the determination accuracy of (3 just like for the other parameters. For the case of xbh = 0.5, we 
obtained the constraint on wbd as wbd > 4854 on average. Comparing this to the one obtained when we set xbh = 0, 
wbd > 4862, it can be seen that the offset of spin parameter affects the constraint on wbd only a little. 

Massive Graviton theories. - 

We fix the masses of the binary constituents to IO^A/q and IO^Mq. We assume the spins of the IO^Mq black holes 
to be for simplicity. As in the case of Brans-Dicke theory, we take the spins of heavier black holes to be x = 0.5. 
Under this assumption, we can apply the simple precession approximation. The results are shown in the second half 
of Table V. The third row in this table represents the results without taking eccentricity le into parameters and the 
fourth row shows the ones including /g. The corresponding histograms are shown in Figs. 13, 14 and 15. The (purple) 
dotted thick line represents the estimation without taking 1^ into parameters and the (red) solid thick one is the 
estimation with /g. In Fig. 13, it is shown that the lower bounds on Xg increase by one order of magnitude when we 
include precessional effect. Compared to the result of [25], which states that the detection of (10^ + 10^)Mq BH/BH 
inspiral gravitational waves with LISA leads to the constraint \g > 1.33 x lO^^cm on average, our result shows that 
the inclusion of a, le and the precessional effect in total enhances the constraint by a factor of 2.3. In this case, since 
the parameter estimation only depends weakly on a and le, the important point is to include the precession which 
makes the constraint stronger. Figure 14 shows that the inclusion of precession considerably enhances the accuracy 
of A/x//x and A/3 by 2 or 3 orders of magnitude. In Fig. 15, one can see again that the tails on ADl/Dl disappear 
due to the change in L. Histograms in Figs. 13 and 14 have tails on the lower accuracy side. These correspond to 
binaries with k « 1, where the effects of precession almost reduce to zero. 

The reason why the effect of spin precession is larger for the massive graviton case than the Brans-Dicke one is 
because the additional information other than the phase is crucial for the parameter estimation in the massive 
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graviton case. The phase of restricted 2PN waveform ^'(/) contains five binary parameters but it has only four PN 
terms. This makes these parameters degenerate. This degeneracy is solved to different levels for different binaries 
when we include precession, which makes the dispersion in the histograms broader. On the other hand, in the Brans- 
Dicke case, there are 5 PN terms in the 2PN phase *(/). Therefore, the degeneracies between parameters are not so 
strong and the effect of precession is relatively weak. 



In this paper, we extended the previous analysis by Berti et al. [25] to see how the inclusion of the spin-spin 
coupling (T, the eccentricity eo and the precessional effect affects the binary parameter estimation by means of LISA 
in the context of alternative theories of gravity such as the Brans-Dickc and massive graviton theories. For the Brans- 
Dicke case, we assumed that we detect NS/BH inspiral gravitational waves with SNR=v'200 because the binaries 
composed of different types of compact objects enhance the gravitational dipole radiation. On the other hand, we 
thought of BH/BH inspiral gravitational waves at 3Gpc for the massive graviton case because the remarkable difference 
between Wph and c appears in lower frequency GWs. 

First, we performed the analysis using the pattern-averaged waveform. For the error estimation in Brans-Dicke 
theory, the inclusion of both a and le into parameters reduces the determination accuracy by an order of magnitude. In 
particular, including 7e affects the estimation more than including just a. However, if we impose the prior information 
A/e > 0, the constraint becomes stronger and would be as good as the one without including Jg into parameters. For 
the analysis in the massive graviton theories, the inclusion of these parameters only changes the results by a factor of 
a few. In this case, the inclusion of a affects more than the inclusion of le- Also in the massive gravity case, imposing 
the prior distribution A/g would not affect the constraint on \g. 

Next, we performed the Monte Carlo simulations including precession. We set the dimensionless spin parameter x of 
the lighter body of binaries to and that of the heavier body to 0.5, which corresponds to setting the Kerr parameter 
to the half of mass in the black hole case. We found that in Brans-Dicke case, the results are not so much affected 
by taking precession into account. For a NS/BH binary of (1.4 -|- 1OOO)M0 with SNR=-\/200, estimation with taking 
a and precession into account can Iciad to a constraint Wbd > 6944 on average (/g is not important for the constraint 
on cjbd because of the prior information for /g). However, when we consider the event rate, the dctectability of 
(1.4-1- 1OOO)M0 binaries of SNR=V200 is very low. Therefore only a lucky detection can constrain wbd- For example, 
referring to Brown et al. [54], the event rate for EMRI is given by 0.7 x (SOOMq/M) x 10"i"Mpc"^yr"^ The distance 
of (1.4 + 1000)Mq binaries of SNR=\/200 approximately corresponds to 40Mpc for LISA. Therefore the event rate 
is 7.6 X lO^'^yr"^ However, for DECIGO and BBO, the event rate of NS/stellar mass BH binaries with SNR=V200 
is about 10^/yr (see [55] and references therein). Therefore, these binaries are thought to be the definite sources for 
them. By using these detectors, we will obtain stronger constraint [56] , because (i) the number of gravitational cycles 



is larger, (ii) the velocity at 1 yr before coalescence is slower and (iii) the width of effective frequency range of 
observation is larger. We can further improve our constraints by using statistical opportunities of large event rates. 
We show these results in a separate paper [56] . 

In the case of massive graviton theory, inclusion of precession has more remarkable effect. The constraint on Ag 
now becomes an order of magnitude stronger. For a BH/BH binary of (10^ 4- 1O^)M0 at 3Gpc, estimation with taking 
cr, le and precession into account can constrain the graviton Compton wavelength as Xg > 3.06 x lO^^cm on average. 
This is 2.3 times stronger than the result obtained in Ref. [25]. 

Our results are consistent with previous results, although our error estimates might be underestimates since our 
analysis does not include systematic errors [38]. 

In Ref. [57], authors have derived binary waveforms in the frequency domain including higher order eccentricity 
terms up to Gq. We included these terms and calculated the constraints on both wbd and Xg with pattern- averaged 
analyses. We found that with our choice of fiducial eccentricities eo = 0.01 at 1 yr before coalescence, the inclusion 
of higher order effects do not affect the constraints. These higher order terms affect the constraints on wbd only 
when eo is larger than 0.1. It seems that this rather high eccentricity is not realised at 1 yr before coalescence from 
the discussion in Sec. VII A and Appendix D. The effect of these higher order terms is much weaker for the massive 
gravity case as the correlation between eccentricity and Xg is weaker compared to the Brans-Dicke case. 

In our calculation, we neglected the spin of one of the binaries. To make our analysis more general, we need to 
take the spins of both binary stars into account. In that case, because the simple precession approximation is no 
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longer valid anymore, we need to solve the precession equations as Stavridis and Will did [39]. We can also improve 
our analysis by introducing higher harmonics to the waveform. Recently, Arun and Will [58] estimated the possible 
constraint on Xg using the higher harmonics waveform, assuming that they detect the SMBH/BH inspiral GWs with 
LISA, advanced LIGO and ET. Since they did not take spins into account, we should try to include both precession 
and higher harmonics into analysis together, and perform the Monte Carlo simulations in alternative theories of 
gravity such as the Brans-Dicke and the massive graviton theories. 
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Appendix A: DETECTED SIGNALS 

The interferometer having three arms corresponds to having two individiial detectors. Therefore, it is possible to 
measure both polarisations with one detector. We first focus on the detector I which consists of the arms 1 and 2. 
This detector measures 

where SLi(t) and 57^2 (i) ^r^' the differences in length in arms 1 and 2 as caused by passing gravitational waves. L is 
the length of the arms where gravitational waves are not present. The factor •\/3/2 comes from the 60° opening angle 
of adjacent arms. Here, we introduce two principal axes for the wave; p = N x L/\N x L\ and q = px N, where L 
is the unit vector parallel to the orbital angular momentum and TV is the unit vector pointing toward the centre of 
mass of the binary. Then, the two polarisations become exactly 7r/2 out of phase and the waveform becomes 

hab(t) = A+H+ COS (f>{t) + A,, H^^ sin ^{t), (A2) 
where ff^^ and H^f^ are the polarisation basis tensors, defined as 

^ab = Po.Pb - QaQb, H^b = P«-1b + laPb- (A3) 

Prom Eq. (Al)-Eq. (A3), the detector output hi{t) becomes 

hit) = ^A+F+{0sAs,i^s)cos^{t) + ^AxFi^(0s,'/>s,V's)sin0(O. (A4) 

F^{9s, (psi^s) and F^ {9s, <^Sj "^s) are the detector beam-pattern coefficients and when the detector is an interferom- 
eter, they are given by 

F+{9s,<ks,iJs) = i(l + cos2 0s)cos(20s)cos(2V's)-cos(^s)sin(2</)s)sin(2?As), (A5) 
Fi''{9s,(t>s,^s) = ^(l + cos2es)cos(2<^s)sin(2V's) + cos(0s)sin(2<^s)cos(2V's). (A6) 
(^Sj (t's) represents the direction of the source in the detector frame and Vs is the polarisation angle defined as 

qz L- z-{L- N){z- N) 

tanti-s = ^-^ = ^-^^-^ — -. (A7) 

N-{Lxz) 
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Also, wc can think of another detector consisting of arms 2 and 3. We cah this detector IF and the signal of this 
detector can be written as hw = {6L2{t) — 6L^{t))/ L. However, since hi and hiji have some correlations, they are not 
independent detectors. We combine detectors I and IF to construct detector II which is uncorrelated with detector I. 
The signal of detector II is 



2 {hxy + hyx) 



(A8) 



This detector II corresponds to an interferometer that is rotated by 45° with respect to detector I. Thus the beam- 
pattern coefficients for the detector II are 



fl^(^s,<^s,V's) 



Fj+(^s,0s-7r/4,Vs), 
j;^(^s,'/'s-7r/4,^s). 



(A9) 
(AlO) 



Reexpressing the waveforms measured by each detector in terms of an amplitude and phase, they become Eq. (19). 

When we perform parameter estimation, we take the direction of the source {Os,(t>s) and the direction of the orbital 
angular momentum (^l, both measured in the solar barycentric frame, as binary parameters. Therefore we need 
to express the waveforms (especially L-N and the beam-pattern functions F+ and which appear in Eqs. (20)-(22)) 
in terms of Ss,4>s, and 4>l- Os{t) and ^s{t) are expressed as 



1 - \/3 
cos^s(t) = 2 ^°^^s - ^ sin6is cos[(f){t) - (/>s], 

TT _i /yScos^s +sin^scos[(^(t) - (^~s] 

Mt) = ai + — + tan „ . . ri... — ^1 

12 \ 2sm6'ssm[0(f) — 0sj 



(All) 
(A12) 



For the polarisation angle Vs (see Eq. (A7)), first z - N = cos^s- Next when we neglect the spin precessional effects, 
L is constant and L ■ z, L ■ N, and N ■ {L x z) are given as 



L-z 
L-N 

N-{Lx z) 



1 - 

- cos 9l Y cos[(^(t) - 

COS cos + sin 9^ sin cos(^l — "^s ) j 
i sin Ol sin sin(<^L - <^s ) 
\/3 

— — COS ^{t) (cos sin sin (ps — cos sin sin ) 
\/3 

— — sin ^{t) (cos sin Oh cos ^i, — cos sin 6s cos ) • 



(A13) 
(A14) 



(A15) 



Appendix B: THE ANALYTIC FORMULAE FOR THE ORBITAL ANGULAR MOMENTUM L UNDER 

SIMPLE PRECESSION 

The precession equations for circular orbit binaries are [34] 



L = ^ 



Si 

S2 



4m 1 



2mi 



3m2 ^ , 4m2 + 3mi 
iJi H -z iJ2 



2m2 



X L 



■1\[{S2 ■ L)S, + (Si • L)S2\ X i - (^\ L, 

2 r"* 5 r \ r 



1 



4rni + 3m2 

2toi 
4to2 + 3toi 

2m2 



X Si 



1 



X S2 + 



S2 

Si 



{S2 ■ L)L 
{Si ■ L)L 



xSi, 

X S2. 



(Bl) 
(B2) 
(B3) 
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The first term of each equation represents the spin-orbit interactions (1.5PN) and the second term represents the 
spin-spin interactions (2PN). The last term of Eq. (Bl) is the angular momentum loss due to the radiation reaction. 
This changes the total angular momentum J = L + Si + S2 as 



5 r \ r J 

In this paper, we assume that one of the spins of the binary constituents is negligible (i.e. Si ~ 0). Then, 
there do not exist spin-spin interactions. We also assume that the orbital angular momentum L is neither parallel nor 
antiparallel to the total spin angular momentum S{= Si + 82)- Then, the precession equations are simplified and L is 
obtained analytically up to some approximate orders. This is the so-called simple precession approximation [34] . This 
also holds when the masses of the binary constituents are equal (mi 7712) and spin-spin interactions are negligible, 
instead of S'l ^ 0. Under this simple precession approximation, the precession equations become Eqs. (30)- (33). 

Next, we define a quantity j{t) as 



7(0 = (B5) 

Then, the magnitude and the direction of the total angular momentum J can be expressed in terms of k, 7(t), L{t), 
L and S, 



J = LVTT2^^7 + 7^, (B6) 
J = , ^ + ^^ . (B7) 



From Eqs. (32), (33) and (B7), the precession equation of J can be derived as 



y _ j[S{l + Kj)-L{K + 'y)] 

- (1 + 2K7 + 72)3/2 • 

Prom Eqs. (32) and (33), it can be seen that the vectors L and S precess around J with the angular velocity Op 
given as Eq. (34). In general, the precessing time scale f2~^ is shorter than the inspiral time scale L/\L\. Therefore 

from J = LL, J changes in magnitude but the direction is almost constant. Actually, if J is much smaller than L 
(as this can happen when L and S are antialigned with almost the same magnitudes), J can change significantly in 
one processional period. Therefore, we introduce the following small parameter. 

Then, J precesses around the fixed direction Jq with 

Jo = J - eJ X L (BIO) 
up to 0{e^). To the same order, the precession equation for L becomes 



L = fipJo X L + eilp{Jo X L)x L. (BH) 

The solution of this equation can be obtained geometrically [34]. We take the baryccntric Cartesian frame (x,y,z) 
which is tied to the ecliptic and centred in the solar system barycentre. Let the barycentre of the binary point in the 
((^s, ^s)) direction and let Jo point in the (^j, ^j) direction. We denote Al as the opening angle of the cone on which 
L processes (i.e. the angle between L and Jq; see Fig. 17). This can be regarded as the angle between L and J apart 
from the errors of order and is given by 
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FIG. 17: The unit orbital angular momentum vector L precesses around a fixed vector Jq. The opening angle of the precession 
cone is given by Al and the precession angle is denoted as a. 



COS At 



L J 



1 + K7 



sin At, — — 



\L\ _ 7VI - 



+ 2k7 + 72 



Then, L can be expressed as 



(B12) 
(B13) 



f i Z-Jocos9,^ ~ ^ sinALsina 

L = Jq cos Al H ; — ^ sin Al cos a + [Jo x z) 



sin 6'j 

where a is the precession angle defined as the solution of 



sin 6 J 



(B14) 



dt - 

We assume that a = is when L ■ z \s maximum (see Fig. 17). By solving the above equation, we get 



(B15) 



5 1/ 3 TOi 

a =Q.r T. T 1 H 

96 y?M^ V 4 ma 



1(CiLf - 3kS{L + KS)gL - 3kS'^(1 - K2)arcsinh 



where ac is a quantity which characterises a at t = tc and Q is defined as 



f L + kS 



(B16) 



g = V1 + 2k7 + 7^. 



(B17) 



From Eq. (B14), the quantities {L ■ N),{L ■ z) and [N ■ {L x 2)], which are needed to compute the polarisation 
angle ips (t) in the beam-pattern coefficients and , are expressed as [28] 
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L-z 



I i X , l-2(Jo-i)cos6ij . 

= 1 Jo • 2;) cos Al H ^ smALCOsa 

2 sin O3 



+— J- Sin Al Sin a, (B18) 

sin 9 J 



r Kr (i ^n \ , cos - (io ■ J^j cos . 

L-N = (Jo ■ N) cosXl -\ = smALCOsa 

sin^j 

(Jo X I) ■ AT . . 
: — = sin Al sin a. 



H = sin Al sin a, (B19) 

sin 0] 

<r ,t- -.^ tCt , i -.^ ^ iV • (I X i) - TV • (Jo X i) cos^j . , 

N -{Lx z) = N ■ {JqX z) cos Al + ^ sm Al cosa 

sin 6.1 



sin 6 J 

N ■ (Jo X z) X £ , X 

+ ^ " - ' sin Al sin a, (B20) 

sin^j 



where 



1 - 

Jo-i = -cos6lj —sm9jcos[(l>{t) - (j)]], (B21) 

Jq • iV = cos^j cos^s + sin^j sin^s cos(^j — ^s), (B22) 
JV • ( Jo X z) = ^ sin 6 J sin sin(^j - ^s) 

— — cos (/)(t) (cos 9 J sin sin 0s — cos sin 9j sin 0j ) 
\/3 

— — sin 0(t) (cos sin cos 0j — cos 9j sin cos ^s), (B23) 

( Jo X z) • iV = sin^s sin^j sin(0j - 0s), (B24) 
/3 - - 

iV-(lxz) = ^sin^ssin(0(f)-0s), (B25) 

(Jo X i) • z = — sin^j sin(0(t) - 0j), (B26) 

iV-(JoXz)xz = -^sin^j[\/3cos^scos{0(t) - 0j} + sin^scos(0j - 0s)- (B27) 

When spins are zero, sinAL = 0, = ^l,<^j = 4>l- Then, Eqs. (B18), (B19) and (B20) each reduces to 
Eqs. (A13), (A14) and (A15), respectively. 

L ■ u and iV • (X x w) in Eq. (46) are given as follows: 

L u = -TV • (L X Jo) (B28) 

- ,± ,,sinAL , „ '..sinAi . 

= TV • (Jo X z) —cosa — (cosys — cosc'jf Jo • TV)) — sina, 

sin 9 J sin 

N-{Lxu) = cosAL-(i-TV)-(Jo-TV), (B29) 
where TV • (Jq x z) and Jo • TV are given as Eqs. (B23) and (B22), respectively. 

Appendix C: THE INVERSION OF THE FISHER MATRIX 

If the ratio (we denote this by R) of the smallest eigenvalue to the largest one in the Fisher matrix T approaches 
the machine's floating-point precision, the inverse of T will not be performed correctly. This problem can be avoided 
by the following technique. 
First, we define 
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diag 



(CI) 



Next, we obtain the normalised Fisher matrix F' as fonows: 



r' = TFT^^ = TFT. (C2) 

Then, all the diagonal components of F' equal to 1. After that, we take the inverse of F'. From F'~^ = x~^F~^T~^, 
we obtain the inverse of our original Fisher matrix F by multiplying T from both sides of F'~^, 



p-i ^ xF'-^T. (C3) 

Even if the ratio R is smaller than the machine's floating-point precision, the ratio for the normalised Fisher matrix 
F' can be larger than the floating-point precision so that the inversion can be performed correctly even in the double 
precision computation. 

To check that our inversion is correctly performed, we followed a similar procedure described in the Appendix of 
Berti et al. [25]. We simply multiply the inversed Fisher matrix by the original one so that we obtain a numerical 
"identity matrix" 7?"™. Then, we define the following small quantity 



(C4) 

With our inversion, we found that ejnv easily satisfies the criteria, £inv < 10~^ (for NS/BH binaries) and Emv < 10~^ 
(for BH/BH binaries), which are used in Ref. [25]. 



max 

i.i 



rnum- rnum 
ij ji 



Sij 



Appendix D: The evolution of ECCENTRICITY 

In this section, we explain how to calculate eccentricity e at a given frequency /, with some initial eccentricity 
e = e, at a frequency / = /«. According to Peters [27], the evolution equations for the semimajor axis a and the 
eccentricity e of the orbit are given by 



da 

de 
di 

Dividing Eq. (Dl) by Eq. (D2), we get 



64 mim2M 
"Ta3(l-(;2)7/2 

304 niiDi-zM 
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304^ 



(Dl) 
(D2) 



Then, we integrate this equation to give 



da 
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/ \ 12/19 

(i) 
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121 2-1 
304^ 
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.1 + 


121 p2 
304 "^i -1 



870 
2299 



(D4) 



where the subscript i denotes the value at the initial time. When e/ej <C 1, Eq. (D4) becomes 



(D5) 



Next, we calculate the time to coalescence Tc{ai, e^). This can be derived from substituting Eq. (D4) into Eq. (D2), 
taking the reciprocal and performing the integration with e from to e^. 
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J^. de 304mim2M ' (i + i2ie2)5299 Jo (1 - e2)3/2 

In the limit of ^ 0, the equation above reduces to 



Then, Eq. (D6) can be expressed as 



T (a- e.) = (aOe-^^/^^i^^ e-/-ii±ii^rfe (D8) 



(■*■ ' 304 

When ~ 1, this equation can be written as 



I^c(ai,ei)«^r,(aO(l-e2)V'. (D9) 



From these equations, we can calculate the eccentricity e at frequency / for the binary whose initial eccentricity is 

5 

256"' 



1 and coalesces in time t{f) = ■^M{'KMf) (see Eq. (16)). From Eqs. (D7) and (D9), we eliminate Tc{ai) 



to get 



5 (l-e2)V2 



(DIO) 



where Tc(ai,ei) is given by Tc(ai,ei) = t(fi). On the other hand, when e <C 1, we can replace the semimajor axis 
a{f) at frequency / by the orbital separation r{t{f)) as [40] 



«(/) « rm) = i^-^^,Mt{f)j . (Dll) 

Then, we obtain our value e by substituting Eqs. (DIO) and (Dll) into Eq. (D5). For (1.4+103)Mq NS/BH binary 
with initial eccentricity Cj = 0.998 at /j = 2 x 10~^Hz, the eccentricity e at frequency / = /lyr becomes e = 0.026. 



Appendix E: THE ANALYTIC FORMS OF ^ USED IN THE PATTERN-AVERAGED ESTIMATES 



In Sec. VII A I, we calculate the parameter determination errors with pattern-averaged estimates. There are 

9 parameters in total: In Al, In 77, fc. oi'c. -Dl, /3, cr, gq, and ui or Bg. The waveform is given by Eq. (29), where the 
amplitude A and the phase ^'(/) are given by Eqs. (27) and (28), respectively. The derivative of this waveform with 
respect to each parameter is taken analytically as [25] 
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where 



dlnM 128 

+a4X + b4X^^'^ + C4x'^)h, (El) 



dlnrj 96' 

+a5X + hx^/"^ +C5x'^)h, (E2) 
= -h, (E3) 



dh 



8t. = '^""^^^ (^^^ 

i = (E7) 
9/i 5i 



<S2(7rA^/)-5/3a;-^/i, (E9) 

(9/) 

— = -i{TrMf)-^h, (ElO) 



aw 3584 
dh 



4 /743 11 \ 128^ 2/5 

64 = ^(/3-47r), (E12) 

5 

/ 3058673 5429 617 , \ . ^ 

= ^1016064 + 1008'^+ 144^ -^j' ^^^^^ 

^ _?!^7^, (E14) 

1462 ^ ' 

ki = -^w, (E15) 
743 33 

= 168 - T^' (^''^ 

65 = ^(^-47r), (E17) 

5 

„/ 3058673 5429 617 , \ 

= '\immM-m^2^-^^-n^ ^^^^^ 

.^5 = ^/e, (E19) 

^ 5848 ^ ' 



34 



[1] Supernova Search Team CoUab. (A. G. Riess et al), Astrophys. J. 607, 665 (2004). 

[2] E. J. Copeland, M. Sami and S. Tsujikawa, Int. J. Mod. Phys. D15, 1753 (2006). 

[3] Y. Fujii and K. Maeda, The Scalar-Tensor Theory of Gravitation, Cambridge University Press (2007). 

[4] P. J. Steinhardt and F. S. Accetta, Phys. Rev. Lett. 64, 2740 (1990). 

[5] C. Brans and R. H. Dicke, Phys. Rev. 124, 925 (1961). 

[6] B. Bertotti, L. less and P. Tortora, Nature 425, 374 (2003). 

[7] V. A. Rubakov and P. G. Tinyakov, Phys. Usp. 51, 759 (2008). 

[8] M. Fierz and W. PauU, Proc. Roy. Soc. Loud. A173, 211 (1939). 

[9] H. van Dam and M. J. G. Vcltman, Nucl. Phys. B22, 397 (1970). 
[10] V. I. Zakharov, JETP Lett. 12, 312 (1970). 
[11] A. L Vainshtein, Phys. Lett. B39, 393 (1972). 

[12] G. R. Dvali, G. Gabadadzo and M. Porrati, Phys. Lett. B485, 208 (2000). 

[13] A. Nicolis and R. Rattazzi, JHEP 0406, 059 (2004). 

[14] V. A. Rubakov, hcp-th/0407104. 

[15] S. L. Dubovsky, JHEP 0410, 076 (2004). 

[16] C. Talmadge, J. P. Berthias, R. W. HelUngs and E. M. Standish, Phys. Rev. Lett. 61, 1159 (1988). 

[17] K. Danzmann, Class. Quant. Grav. 14, 1399 (1997). 

[18] LISA web page: http://lisa.nasa.gov/ 

[19] D. M. Eardlcy, Astrophys. J. 196, L59 (1975). 

[20] C. M. Will, Astrophys. ,J. 214, 826 (1977). 

[21] C. M. Will and H. W. Zaglaucr, Astrophys. J. 346, 366 (1989). 

[22] C. M. Will, Phys. Rev. D50, 6058 (1994). 

[23] P. D. Scharre and C. M. Will, Phys. Rev. D65, 042002 (2002). 

[24] C. M. Will, and N. Yuncs, Class. Quant. Grav. 21, 4367 (2004). 

[25] E. Bcrti, A. Buonamio and C. M. Will, Phys. Rev. D71, 084025 (2005). 

[26] C. M. Will, Phys. Rev. D57, 2061 (1998). 

[27] P. C. Peters, Phys. Rev. 136, B1224 (1964). 

[28] A. Vecchio, Phys. Rev. D70, 042001(2004). 

[29] R. N. Lang and S. A. Hughes, Phys. Rev. D74,122001 (2006). 

[30] M. V. van der Sluys, C. Roever, A. Stroeer, N. Christensen, V. Kalogera, R. Meyer and A. Vecchio, Astrophys. J. 688, 
L61 (2008). 

[31] M. V. van dcr Sluys, V. Raymond, I. Mandcl, C. Rover, N. Christensen, V. Kalogera, R. Meyer and A. Vecchio, Class. 

Quant. Grav.25, 184011 (2008). 
[32] V. Raymond, M. V. van der Sluys, I. Mandel, V. Kalogera, C. Rover and N. Christensen, arXiv:0812.4302 [gr-qc]. 
[33] M. V.van der Sluys, I. Mandel, V. Raymond, V. Kalogera, C. Roever and N. Christensen, arXiv:0905.1323 [gr-qc]. 
[34] T. A. Apostolatos, C. Cutler, G. J. Sussman, and K. S. Thorne, Phys. Rev. D49, 6274 (1994). 
[35] L. Blanchet, T. Damour, B. R. Iyer, C. M. Will and A. G. Wiseman, Phys. Rev. Lett. 74, 3515 (1995). 
[36] L. Barack and C. Cutler, Phys. Rev. D69, 082005 (2004). 

[37] Just after we have completed our analysis, there appeared a paper by Stavridis and Will [39], in which they estimated the 
possible bound on graviton mass from SMBH/BH inspiral gravitational waves using LISA. Following Lang and Hughes [29], 
they randomly distributed the spins of both binary constituents and solved the precession equations numerically. They took 
spin-spin coupling a into account but they did not include eccentricity into parameters. When neglecting the precessional 
spin-spin coupling effect, they obtained the probability distribution of \g which peaks around Ag = 4 x lO^^cm a BH/BH 
binary of (lO" + 1Q^)Mq at 3Gpc. To compare our calculation with their result, we performed another Monte Carlo 
simulation of (10^ + 1.1 x 1O®)M0 binaries. We set the spin of one of the black holes to and we randomly distribute the 
spin of the other black hole from x = to 2. (x is the dimensionless spin parameter defined in Sec. II. Although x > 1 is 
unphysical, we extended it up to x = 2 just to compare our result with the one in Ref. [39].) We included the spin-spin 
coupling a into binary parameters but did not take eccentricity into account. We obtained the probability distribution 
which has a peak around \g = 3.7 x 10^^ cm. Although we cannot directly compare these two results, it seems that our 
result is consistent with the one obtained by Stavridis and Will [39]. 

[38] C. Cutler and M. Vallisneri, Phys. Rev. D76, 104018 (2007). 

[39] A. Stavridis, C. M. Will, arXiv:0906.3602 [gr-qc]. 

[40] C. Cutler and E. E. Flanagan, Phys. Rev. D49, 2658 (1994). 

[41] C. Cutler and J. Harms, Phys. Rev. D73, 042001 (2006). 

[42] A. Krolak, K. D. Kokkotas and G. Schaefer, Phys. Rev. D52, 2089, (1995). 

[43] C. Cutler, Phys. Rev. D57, 7089 (1998). 

[44] L. S. Finn, Phys. Rev. D46, 5236 (1992). 

[45] A. J. Farmer and E. S. Phinney, Mon. Not. Roy. Astron. Soc. 346, 1197 (2003). 

[46] G. Nelemans, L. R. Yungelson, and S. F. Portegies Zwart, Astron and Astrophys. 375, 890 (2001). 

[47] S. A. Hughes, Mon. Not. Roy. Astron. Soc. 331, 805 (2002). 

[48] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes in Fortran, Cambridge University 



35 



Press (1992). 

[49] C. Hopman and T. Alexander, Astrophys. J. 629, 362 (2005). 
[50] P. J. Armitage and P. Natarajan, Astrophys. J. 634, 921 (2005). 

[51] I. Berentzen, M. Preto, P. Berczik, D. Merritt and R. Spurzem, Astrophys. J. 695, 455 (2009). 
[52] M. Vallisneri, Phys. Rev. D77, 042001 (2008). 

[53] B. Kocsis, Z. Haiman, K. Menou and Z. Prei, Phys. Rev. D76, 022003 (2007). 

[54] D. A. Brown, J. Brink, H. Fang, J. R. Gair, C. Li, G. Lovelace, I. Mandel and K. S. Thorne, Phys. Rev. Lett. 99, 201102 

(2007). 

[55] M. Shibata, K. Kyiitoku, T. Yamarnoto and K. Taniguchi, Phys. Rev. D79, 044030 (2009). 
[56] K. Yagi, and T. Tanaka (in preparation). 

[57] N. Yunes, K. G. Arun, E. Berti, and C. M. Will, Phys. Rev. D80, 084001 (2009). 
[58] K. G. Arun and C. M. Will, arXiv:0904.1190 [gr-qc]. 



